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ABSTRACT 



We address the question of to what accuracy remote sens- 
ing images of the surface of planets can be matched, 
so that the possible displacement of features on the sur- 
face can be accurately measured. This is relevant in the 
context of the libration experiment aboard the European 
Space Agency's BepiColombo mission to Mercury. We 
focus here only on the algorithmic aspects of the problem, 
and disregard all other sources of error (spacecraft posi- 
tion, calibration uncertainties, etc.) that would have to be 
taken into account. We conclude that for a wide range of 
illumination conditions, translations between images can 
be recovered to about one tenth of a pixel r.m.s. 
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1. INTRODUCTION 



One of the goals of the European Space Agency's Bepi- 
Colombo mission to Mercury is the measurement of the 
amplitude of the libration of Mercury. In order to do this 
images of the same surface areas will be taken at differ- 
ent times during the libration cycle and compared. When 
all other effects — spacecraft position. Mercury's rotation, 
spacecraft attitude, etc. — are taken into account, any re- 
maining discrepancy between the positions of features on 
the surface must be due to the libration of the crust of the 
planet. 

Here we address the question of to what accuracy images 
can be matched, and we focus only on the algorithmic 
aspects of the problem, disregarding all other sources of 
error that would have to be taken into account to solve 
the scientific problem. We shall show that by using a 
shape-based matching algorithm images taken under a 
wide range of illumination conditions can be matched to 
one tenth of a pixel root-mean-square. Based on this we 
conclude that the accuracy of the pattern matching algo- 
rithm is not the limiting factor in the ultimate accuracy 
that can be achieved by the libration experiment on Bepi- 
Colombo. 



2. PATTERN MATCHING 



The pattern matching algorithms to be used in this study 
will have to deal with images that may appear to be dras- 
tically different from one another, still they refer to the 
same region. Consider for example the images shown in 
figure[n To the human eye it is clear that the images refer 
to the same region, but any algorithm that relied on the 
presence of identical features in the images would have 
great difficulty concluding that the images are related at 
all. 

What is clear by visual inspection is that a number of 
edges — sharp changes in the level of illumination be- 
tween contiguous pixels — are common between images. 
These edges appear where sharp changes in the altimetric 
profile occur. It is also clear that not all edges appear in 
all images, owing to the complex interplay between the 
position of the Sun, and the orientation and slope of the 
features on the ground. 

Compare for instance images b and d in figure^ the left 
rim of the crater is bright in one image, and dark in the 
other. No similar change is observed on the right rim of 
the crater 

Take now images a and c. Here the left rim of the 
crater appears almost to be the same, but the extent of 
the shadow cast by the right rim is dramatically different. 

The ideal algorithm must be able to identify the edges in 
the two images, must be robust against local, non-linear 
changes in illumination conditions, and it must be able to 
operate by identifying a subset of features that are com- 
mon to the pair of images being compared. Finally, based 
on the common features identified, the algorithm must be 
able to recover a possible translation between the two im- 
ages. 

Algorithms that try to minimize the difference between 
the two, possibly scaled, images are clearly not going 
to be suited for the task, unless the images to be com- 
pared are taken under very similar illumination condi- 
tions. While this is possible, it would be a very strong 
constraint on the operations of a mission. 
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Figure 1. A digital elevation model of Olympus Mons 
viewed by an imaging camera under different illumina- 
tion conditions: (a) The Sun is at 5° elevation; (b) The 
Sun is at 25° elevation; (c) The Sun is at 50° elevation; 
(d) The Sun is at 85° elevation. In all cases the Sun's 
azimuth is 0° (to the right). 

Based on the considerations above, we have chosen to 
make use of the image matching algorithms available in 
the HALCON software library (Ref. 1). This is a com- 
mercial product used in image vision and image recogni- 
tion applications. 

One particular technique available in the HALCON U- 
brary is the so-called shape-based matching (Ref. 
This technique is based on an algorithm that identifies the 
shape of patterns in images, and can be instructed to find 
in a comparison image the shape identified in a reference 
image. 



2.1. Shape-Based Pattern Matching 

The detailed description of the algorithm can be found in 
the HALCON documentation (Ref. 2) and has been sub- 
mitted as part of a European Patent Application (Ref.0). 

The algorithm proceeds through the following steps: 

1. A so-called region of interest is identified in the ref- 
erence image. This is a region of the image where 
edges will be looked for. The region of interest must 
be selected to be fully contained in both images. 
This step is done by hand, based on some a priori 
knowledge, or visual inspection of the images. In 
our case, where the simulated translations amount 
to a few pixels along either or both the X and Y axes 
(see §|5|i, the region of interest is the whole reference 
image, minus a few pixels around the edges. In the 
case of two partially overlapping images of the same 
region one would choose the intersection of the two 
images. 



2. Features are identified in the comparison image with 
an edge detection algorithm. Pixels identified by the 
edge detection algorithm are part of the reference 
pattern. 

3. The edge detection algorithm is run on the compar- 
ison image. This results in a second collection of 
pixels, the search pattern. 

4. The algorithm now overlays the reference pattern on 
the search pattern. The reference pattern is stepped 
over the search pattern in an attempt to maximize the 
number of overlapping pixels. In doing so the algo- 
rithm is allowed to reduce the number of pixels in 
the reference pattern. The maximum fraction of the 
search pattern that can be discarded in the process 
can be set by the user In our application the refer- 
ence and search patterns can differ vastly. Therefore 
we have allowed the algorithm to throw away up to 
70 % of the pixels. In trying to maximize the over- 
lap between the two patterns, the algorithm can be 
instructed to allow for a rotation and a scaling fac- 
tor. 

5. The algorithm reports the recovered translations and 
the fraction of the pixels in the reference pattern that 
was used to find a match. The latter is called the 
score. Within the parameters given by the user, the 
algorithm always chooses the match with the highest 
score. 



2.1.1. The Meaning of the Score 

The HALCON score is the normalized sum of the cross 
product between the vectors describing the position of the 
pixels in the reference pattern and those describing the 
pixels in the search pattern. If the two patterns are iden- 
tical, it is clear that the score will be equal to one. When 
pixels have to be dropped from the reference pattern, the 
score will decrease. 

In the actual algorithm, the sum of the cross products of 
the pixels used in the match is slightly modified to take 
into account the possibility of non-linear changes in the 
illumination conditions, either locally, around certain fea- 
tures, or globally, across the entire image (Ref.U. 

It is tempting to interpret the HALCON score as a qual- 
ity factor for the goodness of the translation parameters 
found. However this would be wrong on two counts. 

First of all, it is clear that often only a subset of the pat- 
tern to be looked for is to be found in the search pattern. 
(Refer back to the examples shown in figure In this 
case the search algorithm must discard some of the pixels 
in the reference pattern in order to find a good matching 
sub-pattern. How many pixels are left in the sub-pattern 
has nothing to do with whether the match is good or not. 

Second, the notion of goodness of match implies that the 
matched pattern can be compared to an expected result. 



or true pattern, or that the algorithm proceeds through the 
optimization of an objective function. But the only mea- 
sure of how well the algorithm has performed, is how 
close the recovered translation is to the values injected in 
the simulation. This means that the accuracy of the al- 
gorithm can only be judged through an extensive set of 
Monte Carlo simulations. Only by repeatedly comparing 
two images in multiple realizations of the same detection 
and matching process, is it possible to gauge the statis- 
tical errors in the results, and therefore establish to what 
accuracy and under what conditions the algorithm can be 
effectively employed. 



3. APPROACH 



The following steps were identified. 

1 . Render a digital elevation model of the surface of a 
planet, by choosing the position of the Sun and of 
the spacecraft. We have used povray (Ref. 4) for 
this task. 

2. Create two images, the reference image and the 
comparison image. The latter is possibly translated 
along one or both of the image axes. 

3 . Convert the rendered images to instrument count im- 
ages of a given signal-to-noise ratio. 

4. Recover translation parameters between the two im- 
ages using a shape-based pattern matching algo- 
rithm. 

5. Study the accuracy with which the parameters are 
recovered, and derive information on the range of il- 
lumination conditions for which the parameters can 
be successfully recovered. 



4. DIGITAL ELEVATION MODELS 



Four digital elevation models have been used in this 
study. These are shown in figure|2j- 



5. SIMULATION RUNS 

After some preliminary simulations used to determine a 
useful sampling scheme of the parameter space, the bulk 
of the simulations were carried out with the following pa- 
rameters values: 

• Translations: 100 meter in X, in Y, and in X and Y. 

• Sun elevation angles: 10°, 30°, 60°, 90°. 

'The Olympus Mons digital elevation model was kindly provided to 
us by the Mars Express Team. 
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Figure 2. The four digital elevation models used in this 
study. From left: The Olympus Mons caldera; a bowl- 
shaped crater close to the Olympus Mons caldera; a syn- 
thetic landscape with several bowl-shaped craters; an- 
other synthetic landscape containing approximately 5000 
craters. (An image of Mercury taken from an height of 
400 km might contain a few thousands of craters with a 
diameter larger than a few tens of meters.) Darker colors 
represent lower elevations. 

• Nominal spacecraft height 1500km^. 

• Sun azimuth angle: several (the same azimuth an- 
gle for reference and comparison images). For one 
model a difference in azimuth of 30° between the 
two images was introduced. 

• Four digital elevation models rendered with a signal- 
to-noise ratio of 50. The signal-to-noise ratio is de- 
termined when the Sun is at the zenith. 



6. RESULTS 



For each digital elevation model used, several thousand 
data points have been calculated. Each data point refers 
to a particular combination of Sun elevation and azimuth 
for the reference and comparison images, and a transla- 
tion along one or both of the image axes. For each combi- 
nation of parameters, the same number of simulation runs 
(ten) was carried out. 

In the following we use and Aj, to indicate the differ- 
ence between the amplitude of the translation recovered 
by the algorithm and the amplitude of the translation used 
in the simulation. Therefore the expectation value of A^; 
and Aj^ is always 0, and the width of their distributions is 
a measure of the statistical error in the reconstruction. 

In table^we give a summary of how successful the algo- 
rithm has been. For each digital elevation model we give 
the number of realizations (all Sun angles and all trans- 
lations), how many times the algorithm failed to return a 
match, and how many times the returned result was more 
than 2 pixels away from the expected result. The latter 
figure has no special meaning, but is meant to give an 
idea of the global behavior of the algorithm. 

One thing is immediately apparent: for the synthetic 
digital elevation model the algorithm always returned a 

^The actual height of the camera above the surface is not important 
for the results of this study, at least as long as the images recorded from 
different heights show the same level of detail. 
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DEM 


Total 


No match 


Aa; > 2 or Aj, > 2 


A 


9878 


2.8% 


2.3% 


B 


7300 


3.2% 


4.2% 


C 


6254 


0% 


12.6% 



Table 1. Global success statistics for the simulation runs. 
For each digital elevation model (DEM) the following 
data are reported: the total number of independent re- 
alizations; the fraction of realizations for which the algo- 
rithm was not able to find any match; the fraction of re- 
alizations for which either \ or |Ajy| were larger than 
two pixels. DEM keys: A = Olympus Mons, B = bowl 
crater, C = synthetic. 

match, but a larger fraction of the returned answers was 
significantly wrong. Because the synthetic model is sig- 
nificantly more regular than the other two — in particular 
the craters are identical but for a scale factor, the algo- 
rithm has an easier job at finding some matching pattern, 
although relatively more often the pattern found is not the 
good one. 

Based on these observations, we present the results for 
the synthetic model separately. However, we will be able 
to show that the same conclusions on the accuracy of the 
algorithm can be reached for all digital elevation models 
by applying the same selection criteria on the illumina- 
tion conditions. 

In the following sections we use the following notation: 
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Figure 3. The effect of the Sun elevation cut on the dis- 
tribution of Ax (top) and Ay (bottom) for the Olympus 
Mons data. 



• Ocut refers to the following selection: Osun > 10° 
and 9sun 7^ 90° , where dsun is the Sun elevation 
angle in the reference or the comparison image. 

• 4>cut refers to the following selection: |0stiri ~ n x 
90° I > 20 for n G 0, 1,2,3,4, where (^5„„ is the 
Sun azimuth — this is the same in both the reference 
and the comparison images. 



6.1. The Olympus Mons Models 

Most of the high-deviation points come from images 
where the Sun elevation is equal to 90°, or images where 
the Sun elevation is lower than 10° (Ocut)- This is shown 
in figure |3] 

In figure @] we plot the data with 9 cut applied versus the 
Sun azimuth. We observe that the mean of A^^ and A^ 
vary with (j>sun in a quasi-periodic fashion. What we 
observe is that the deviation is larger when the Sun az- 
imuth is orthogonal to one of the image axes. Namely, the 
largest deviations for A^; are observed when (psun ~ 90° 
or 270°, whereas the largest deviations for Aj^ are ob- 
served when (j>sun ~ 0° and 180° . The direction defined 
by the Sun azimuth appears to be a preferential direc- 
tion: translations along this direction can be more accu- 
rately recovered, because the features on the terrain create 
sharper shadows along the direction to the Sun. 
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Figure 4. The average A^ and Ay versus the Sun azimuth 
for the Olympus Mons data. The error bar on each point 
represents the root-mean-square. 
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Based on the data in figurel^we can devise a selection cri- 
terion for the Sun azimuth, so that translations along both 
axes can be recovered with comparable accuracy. The 
selection criterion is that the Sun azimuth must be more 
than 20° away from both image axes. The distributions 
of Aa; and Aj, when both 9cut and (pcut are applied are 
shown in figure [S] 

Figure |5l represents the end point of our analysis. We 
observe that the two distributions are centered on 0, and 
have a width of w 0. 1 pixel root-mean-square. The two 
distributions have a tail in the direction of the translation 
applied in the simulations (-100 m along the X axis, and 
H-lOO m along the Y axis). The nature of this slight asym- 
metry in not understood at present. 
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Figure 6. The distribution of iS.^ (top) and /S.y (bottom) 
for the Olympus Mons data for date where the Sun az- 
imuth of the comparison and reference images differ by 
30°. 
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Figure 5. The distribution of o-nd ^yfor the Olympus 
Mons data, once both the Sun cuts are applied. 



6.2. The Synthetic Model 

As already hinted to, the results based on the synthetic 
digital elevation model give a slightly different picture, 
although the main conclusions do not change. 

The dcut criterion is still effective in rejecting data points 
that return a large deviation from the expectation. 

A point of discrepancy with respect to the Olympus Mons 
data is the behavior of the recovered translations as a 
function of Sun azimuth. Figure shows that the effect 
observed for the Olympus Mons models is almost not ob- 
served here. After the 9 cut criterion is applied, any re- 
maining offset is smaller than 0.05 pixel. 



6.1.1. Changing the Sun Azimuth 



and A J,. 



Finally, figure |8] shows the distribution for A 
Again, the translation is recovered with an accuracy of 
w 0. 1 pixel root-mean-square, but the details of the distri- 
butions differ from what was observed before. 



The bulk of the simulation runs was carried out with the 
same Sun azimuth for both the reference and comparison 
images. We however also made a set of simulation runs 
where the azimuth of the Sun in the comparison image 
was 30° away from the azimuth used in the reference im- 
age; only a translation of 100 meters along the X axis was 
applied. The results are shown in figure |6l Even in this 
case the algorithm is able to recover the injected transla- 
tion with an accuracy of « 0. 1 pixel root-mean-square. 



7. CONCLUSIONS 



We have performed a study of the accuracy with which 
a shape-based pattern matching algorithm can identify 
translations between remote sensing images of the same 
planetary features. We have applied the algorithms in a 
Monte Carlo fashion to digital elevation models (both real 
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Figure 7. The average (top) and (bottom) versus 
the Sun azimuth for the Synthetic Model data. The error 
bar on each point represents the root-mean-square. 



and synthetic) in order to investigate the statistical perfor- 
mance of the procedure. 

We find that for a broad range of illumination conditions 
translations between images can be recovered with an ac- 
curacy of 0.1 pixel r.m.s. 

The algorithm performs best for translations along the 
projected direction to the Sun on the image plane. This 
study shows that translations along both image axes at 
the same time can be recovered with the same accuracy 
of 0.1 pixel as long as the projected direction to the Sun 
lies more than k, 20° away from the same image axes. 

Finally, this study demonstrates that the images to be 
compared need not be taken under the very same illumi- 
nation conditions in order to be effectively matched. For a 
given Sun azimuth, any pair of images taken with Sun el- 
evation angles larger than 10° can be used; images taken 
when the Sun is at the zenith must also be avoided. The 
range of useful illumination conditions is further broad- 
ened because this study concludes that differences in Sun 
azimuth of at least 30° do not affect the accuracy of the 
matching algorithm. 

The error contributed by the matching algorithm is but 
one of the several error contributions to be taken into ac- 
count during the analysis of the data pertaining to the 
measurement of the possible libration of the surface of 
Mercury. This study shows that the accuracy of the pat- 
tern matching algorithm is not a limiting factor in the 
ultimate accuracy of the libration experiment aboard the 
BepiColombo mission to Mercury. 
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Figure 8. The distribution of Ax (top) and Ay (bottom) 
for the Synthetic Model data, once both the Sun cuts are 
applied. 



